#include "sweFunctions.hpp"
#include <cmath>
#include <cassert>

namespace hms
{


/** @brief Returns the maximum stable time step for the given velocity, water
 * depth and cell parameters.
 * 
 * Implements Simons (2019) Eq. 3.15 with Courant No. (Co) set to 1, so that
 * the return value can be multiplied with the desired max Courant No:
 * dt = Co * 4 * A / ( P * ( |v| + sqrt(g*d) ) )
 * A is cell area, P is cell circumference, v is velocity, d is water depth */
scalar courantTimeStep(
	scalar area, scalar circumf, const Ref<const Array2s>& v, scalar d
){
	assert( d>=0 );
	return
		/* Courant number, solved for delta t, see Simons (2019) Eq. 3.15 */
		4 * area / circumf / ( v.matrix().norm() + wavespeed(d) );
}


/** @brief Returns the Chezy coefficient calculated from water depth and
 * Manning's n.
 * 
 * Equation: C = d^(1/6) / n (Simons, 2019: Eq. 2.5) */
scalar chezy( scalar d, scalar n ){
	assert( n > 0 && "Mannings n must be strictly positive!" );
	return pow(d, 1./6) / n;
}


scalar froude( scalar vNorm, scalar d ){
	assert( d >= 0 && "Water depth must be positive!" );
	return vNorm/sqrt(phys::g * d);
}


Array1Xs froudes( 
	const Array2Xs& velocities,
	const Array1Xs& waterdepth
){
	assert (velocities.cols() == waterdepth.cols() && "number of entries in velocities must be equal to number of entries in waterdepths!");
	
	Array1Xs fr {
		velocities.colwise().norm() / sqrt(phys::g * waterdepth)
	};

	return fr;
}


} // end namespace hms
